function PV_ldgh(grid,base,lambda,u,q,k,l,outfile,test_name)
% PV_ldgh(grid,base,lambda,u,q,k,l,outfile,test_name)
%
% Write a (lambda,u,q) solution of the LDG-H method in a format
% readable by paraview (.vtk).
%
% To avoid matching the nodes of the hybrid mesh with those of
% the primal and dual variables, the output is divided in two files.
% 1) the file outfile-lambda.vtk contains the hybrid variable,
% represented as a fully discontinuous field;
% 2) the file outfile-uomega.vtk contains the primal and the dual
% variables, structured as discontinuous fields and
% interpolated on the same grid.
%
% See also: interpolate.

 % First the hybrid variable

 % Interpolate lambda to the required degree
 [Xi,Ti,L,cell_type] = el_resample(grid,base, ...
   reshape(lambda,[base.nk,grid.ns]),k,l,'hybrid');
 % add th third component for 2D results
 if(size(Xi,1)==2)
   Xi(3,:,:) = 0;
 end

 % Initialize the output file
 Xpos = find(outfile=='.');
 Xpos = Xpos(end); % take the last dot in the file name
 fid = fopen([outfile(1:Xpos-1),"-lambda",outfile(Xpos:end)],'w');
 
 % Preamble
 line = '<?xml version="1.0"?>';
 fprintf(fid,'%s\n',line);
 line = '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="BigEndian">';
 fprintf(fid,'%s\n',line);
 line = ' <UnstructuredGrid>';
 fprintf(fid,'%s\n',line);

 % Write the grid
 npoints = size(Xi,2)*size(Xi,3); % total number of points
 n_v4e   = size(Xi,2);    % number of vertexes per element
 n_v4c   = size(Ti,1);    % number of vertexes per cell
 n_c4e   = size(Ti,2);    % number of cells per element
 n_cells = n_c4e*grid.ns; % total number of cells
 line = [ '  <Piece NumberOfPoints="' , num2str(npoints,'%1d') , ...
                 '" NumberOfCells="'  , num2str(n_cells,'%1d') , '">'];
 fprintf(fid,'%s\n',line);

 line = '   <Points>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Float32" NumberOfComponents="3" Format="ascii">';
 fprintf(fid,'%s\n',line);
 if(grid.d==2) % add z to the data
   Xi(3,:,:) = 0;
 endif
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e ',Xi(:,i,j));
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Points>';
 fprintf(fid,'%s\n',line);

 line = '   <Cells>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="connectivity" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ns
   i_shift = (i-1)*n_v4e;
   for j=1:n_c4e
     fprintf(fid,' %i',i_shift+Ti(:,j)-1);
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="offsets" Format="ascii">';
 fprintf(fid,'%s\n',line);
 i_shift = 0;
 for i=1:grid.ns
   for j=1:n_c4e
     i_shift = i_shift + n_v4c;
     fprintf(fid,' %i',i_shift);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="types" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ns
   for j=1:n_c4e
     fprintf(fid,' %i',cell_type);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Cells>';
 fprintf(fid,'%s\n',line);

 % Write the data (lambda)
 line = '   <PointData Scalars="scalars">';
 fprintf(fid,'%s\n',line);
 in = 1;
 line = ['    <DataArray type="Float32" Name="','lambda', ...
         '" Format="ascii">'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e\n',L(i,j));
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </PointData>';
 fprintf(fid,'%s\n',line);
 
 % Close the grid
 line = '  </Piece>';
 fprintf(fid,'%s\n',line);

 % Close the file
 line = ' </UnstructuredGrid>';
 fprintf(fid,'%s\n',line);
 line = '</VTKFile>';
 fprintf(fid,'%s\n',line);

 fclose(fid);

 % Now the two remaining variables

 % Initialize the output file
 fid = fopen([outfile(1:Xpos-1),"-uomega",outfile(Xpos:end)],'w');
 
 % Preamble
 line = '<?xml version="1.0"?>';
 fprintf(fid,'%s\n',line);
 line = '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="BigEndian">';
 fprintf(fid,'%s\n',line);
 line = ' <UnstructuredGrid>';
 fprintf(fid,'%s\n',line);

 % Interpolate to the required degree
 [Xi,Ti,U,cell_type] = el_resample(grid,base, ...
   reshape(u,[base.pk,grid.ne]),k,l);
 [Xi,Ti,O,cell_type] = el_resample(grid,base, ...
   reshape(q,[base.mk,grid.ne]),k,l,'flux');
 % add th third component for 2D results
 if(size(Xi,1)==2)
   Xi(3,:,:) = 0;
 end
 if(size(O,1)==2)
   O(3,:,:) = 0;
 end

 % Write the grid
 npoints = size(Xi,2)*size(Xi,3); % total number of points
 n_v4e   = size(Xi,2);    % number of vertexes per element
 n_v4c   = size(Ti,1);    % number of vertexes per cell
 n_c4e   = size(Ti,2);    % number of cells per element
 n_cells = n_c4e*grid.ne; % total number of cells
 line = [ '  <Piece NumberOfPoints="' , num2str(npoints,'%1d') , ...
                 '" NumberOfCells="'  , num2str(n_cells,'%1d') , '">'];
 fprintf(fid,'%s\n',line);

 line = '   <Points>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Float32" NumberOfComponents="3" Format="ascii">';
 fprintf(fid,'%s\n',line);
 if(grid.d==2) % add z to the data
   Xi(3,:,:) = 0;
 endif
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e ',Xi(:,i,j));
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Points>';
 fprintf(fid,'%s\n',line);

 line = '   <Cells>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="connectivity" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ne
   i_shift = (i-1)*n_v4e;
   for j=1:n_c4e
     fprintf(fid,' %i',i_shift+Ti(:,j)-1);
     fprintf(fid,'\n');
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="offsets" Format="ascii">';
 fprintf(fid,'%s\n',line);
 i_shift = 0;
 for i=1:grid.ne
   for j=1:n_c4e
     i_shift = i_shift + n_v4c;
     fprintf(fid,' %i',i_shift);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '    <DataArray type="Int32" Name="types" Format="ascii">';
 fprintf(fid,'%s\n',line);
 for i=1:grid.ne
   for j=1:n_c4e
     fprintf(fid,' %i',cell_type);
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);
 line = '   </Cells>';
 fprintf(fid,'%s\n',line);

 % Write the data (u,q,...)
 line = '   <PointData Scalars="scalars">';
 fprintf(fid,'%s\n',line);

 line = ['    <DataArray type="Float32" Name="','u', ...
         '" Format="ascii">'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e\n',U(i,j));
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);

 %line = ['    <DataArray type="Float32" Name="','us', ...
 %        '" Format="ascii">'];
 %fprintf(fid,'%s\n',line);
 %for j=1:size(Xi,3)
 %  for i=1:size(Xi,2)
 %    fprintf(fid,'%.15e\n',Us(i,j));
 %  end
 %end
 %line = '    </DataArray>';
 %fprintf(fid,'%s\n',line);

 line = ['    <DataArray type="Float32" Name="','Q', ...
         '" NumberOfComponents="3" Format="ascii">'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%.15e\n',O(:,i,j));
   end
 end
 line = '    </DataArray>';
 fprintf(fid,'%s\n',line);

 line = '   </PointData>';
 fprintf(fid,'%s\n',line);
 
 % Close the grid
 line = '  </Piece>';
 fprintf(fid,'%s\n',line);

 % Close the file
 line = ' </UnstructuredGrid>';
 fprintf(fid,'%s\n',line);
 line = '</VTKFile>';
 fprintf(fid,'%s\n',line);

 fclose(fid);

%
% % points
% npoints = size(Xi,2)*size(Xi,3);
% line = ['POINTS ',num2str(npoints),' float'];
% fprintf(fid,'%s\n',line);
% for j=1:size(Xi,3)
%   for i=1:size(Xi,2)
%     fprintf(fid,'%f ',Xi(:,i,j));
%     fprintf(fid,'\n');
%   end
% end
%
% % cells
% line = ['CELLS ',num2str(size(Ti,2)),' ', ...
%                  num2str((size(Ti,1)+1)*prod(size(Ti,2)))];
% fprintf(fid,'\n%s\n',line);
% for i=1:size(Ti,2)
%   fprintf(fid,'%i',size(Ti,1));
%   fprintf(fid,' %i',Ti(:,i)-1);
%   fprintf(fid,'\n');
% end
%
% % cell type
% line = ['CELL_TYPES ',num2str(size(Ti,2))];
% fprintf(fid,'\n%s\n',line);
% fprintf(fid,'%i ',cell_type*ones(1,size(Ti,2)-1));
% fprintf(fid,'%i\n',cell_type); % write the last one followed by \n
%
% % Write the data
% line = ['POINT_DATA ',num2str(npoints)];
% fprintf(fid,'\n%s\n',line);
%
% line = 'SCALARS u float 1';
% fprintf(fid,'%s\n',line);
% line = 'LOOKUP_TABLE default';
% fprintf(fid,'%s\n',line);
% for j=1:size(Xi,3)
%   for i=1:size(Xi,2)
%     fprintf(fid,'%f ',U(i,j));
%     fprintf(fid,'\n');
%   end
% end
% fprintf(fid,'\n');
%
% line = 'VECTORS omega float';
% fprintf(fid,'%s\n',line);
% for j=1:size(Xi,3)
%   for i=1:size(Xi,2)
%     fprintf(fid,'%f ',O(:,i,j));
%     fprintf(fid,'\n');
%   end
% end
% fprintf(fid,'\n');
%
% fclose(fid);
  
return


